# c('#304D30', '#4F6F52', '#739072', '#86A789', '#D2E3C8', '#B6C4B6',
# '#DFD0B8', '#DEAA79', '#AF8F6F', '#6F4E37')
calc_annual_trait_RGR <- function(data, trait, RGR = "RGR") {
result <- data.frame()
for (compete in unique(data$Compete)) {
for (year in unique(data$comp_year)) {
subset_data <- data %>%
filter(Compete == !!compete, comp_year == !!year) %>%
select(Nenv, sp, ring.num, Compete, growth_time, Light, Moisture,
AP, trait, RGR, paste0(trait, "_minmax_scaled")) %>%
na.omit()
if (nrow(subset_data) > 5) {
# Fit the mixed-effects model
model_formula <- as.formula(paste0(RGR, " ~ ", trait, " *(Light + Moisture + AP) + (1 | sp) + (1|Nenv)"))
model <- lmer(model_formula, data = subset_data)
# Extract slope (first fixed effect coefficient)
slope <- fixef(model)[2]
TLslope <- fixef(model)[6]
TMslope <- fixef(model)[7]
TPslope <- fixef(model)[8]
# Calculate marginal R² using `MuMIn` package for lmer models
library(MuMIn)
r2m <- r.squaredGLMM(model)[, "R2m"]
r2c <- r.squaredGLMM(model)[, "R2c"]
sp_slopes <- c() # Initialize a vector to store slopes for each species
# Iterate through each species
for (species in unique(subset_data$sp)) {
# Subset data for the current species
species_data <- subset_data %>%
filter(sp == species) %>%
na.omit()
# Ensure there is enough data to fit the model
if (nrow(species_data) > 5) {
# Define the model formula dynamically for the given trait
# and RGR
model_formula <- as.formula(paste0(RGR, " ~ ", trait))
# Fit the linear mixed-effects model for the current
# species
sp_model <- summary(lm(model_formula, data = species_data))
# Extract the slope for the trait (first fixed-effect
# coefficient)
sp_slope <- sp_model$coefficients[2, 1]
# Append the slope to the slopes vector
sp_slopes <- c(sp_slopes, sp_slope)
}
}
# Calculate the mean and standard deviation of slopes
normalized_slopes <- (sp_slopes - min(sp_slopes))/(max(sp_slopes) -
min(sp_slopes))
mean_slope <- mean(normalized_slopes, na.rm = TRUE)
sd_slope <- sd(normalized_slopes, na.rm = TRUE)
# Calculate the coefficient of variation (CV) of slopes
cv_slope <- sd_slope/mean_slope
# Extract p-value for the main effect of the trait
p_value <- summary(model)$coefficients[2, "Pr(>|t|)"]
com_itv <- CV4(subset_data[, paste0(trait, "_minmax_scaled")])
# sp_itv <- mean(tapply(subset_data[, paste0(trait,
# '_minmax_scaled')], list(subset_data$sp), CV4), na.rm = T)
subset_data$new_com <- "Com"
kernel.trait <- TPDs(traits = subset_data[, trait], species = subset_data$new_com) # Kernel density estimates of trait for individuals
functional.indices.trait <- REND(TPDs = kernel.trait)
sp.kernel.trait <- TPDs(traits = subset_data[, trait], species = subset_data$sp) # Kernel density estimates of trait for individuals
sp.fun.indices <- REND(TPDs = sp.kernel.trait)
sp_itv <- mean(sp.fun.indices$species$FRichness)
ol <- dissim(sp.kernel.trait)
dis <- ol$populations$dissimilarity
inter_dis <- mean(dis[upper.tri(dis) | lower.tri(dis)])
mean_trait <- mean(subset_data[, paste0(trait, "_minmax_scaled")],
na.rm = T)
sp_meantrait <- tapply(subset_data[, paste0(trait, "_minmax_scaled")],
list(subset_data$sp), mean, na.rm = T)
# inter_dis <- CV4(sp_meantrait)
RGR_var <- sd(subset_data$RGR, na.rm = T)
RGR_mean <- mean(subset_data$RGR, na.rm = T)
result <- rbind(result, data.frame(Compete = compete, sample_size = nrow(subset_data),
trait = trait, comp_year = year, estimate = slope, TLestimate = TLslope,
TMestimate = TMslope, TPestimate = TPslope, TYestimate = NA, Marginal_R2 = r2m,
Conditional_R2 = r2c, p_value = p_value, com_ITV = com_itv, sp_ITV = sp_itv,
Fric = functional.indices.trait$species$FRichness, Fdiv = functional.indices.trait$species$FDivergence,
mean_trait = mean_trait, inter_dis = inter_dis, RGR_var = RGR_var,
RGR_mean = RGR_mean, sp_slope_cv = cv_slope))
}
}
}
return(result)
}
calc_trait_RGR <- function(data, trait, RGR = "RGR") {
result <- data.frame()
for (compete in unique(data$Compete)) {
subset_data <- data %>%
filter(Compete == !!compete) %>%
select(Nenv, sp, Compete, growth_time, Light, Moisture, AP, trait, RGR,
paste0(trait, "_minmax_scaled")) %>%
na.omit()
subset_data$growth_time <- scale(as.numeric(subset_data$growth_time))
if (nrow(subset_data) > 2) {
# Fit the mixed-effects model
model_formula <- as.formula(paste0(RGR, " ~ ", trait, " *(Light + Moisture + AP + growth_time) + (1 | sp) +(1|Nenv)"))
model <- lmer(model_formula, data = subset_data)
# Extract slope (first fixed effect coefficient)
slope <- fixef(model)[2]
TLslope <- fixef(model)[7]
TMslope <- fixef(model)[8]
TPslope <- fixef(model)[9]
TYslope <- fixef(model)[9]
# Calculate marginal R² using `MuMIn` package for lmer models
library(MuMIn)
r2m <- r.squaredGLMM(model)[, "R2m"]
r2c <- r.squaredGLMM(model)[, "R2c"]
sp_slopes <- c() # Initialize a vector to store slopes for each species
# Iterate through each species
for (species in unique(subset_data$sp)) {
# Subset data for the current species
species_data <- subset_data %>%
filter(sp == species) %>%
na.omit()
# Ensure there is enough data to fit the model
if (nrow(species_data) > 5) {
# Define the model formula dynamically for the given trait
# and RGR
model_formula <- as.formula(paste0(RGR, " ~ ", trait))
# Fit the linear mixed-effects model for the current species
sp_model <- summary(lm(model_formula, data = species_data))
# Extract the slope for the trait (first fixed-effect
# coefficient)
sp_slope <- sp_model$coefficients[2, 1]
# Append the slope to the slopes vector
sp_slopes <- c(sp_slopes, sp_slope)
}
}
# Calculate the mean and standard deviation of slopes
normalized_slopes <- (sp_slopes - min(sp_slopes))/(max(sp_slopes) - min(sp_slopes))
mean_slope <- mean(normalized_slopes, na.rm = TRUE)
sd_slope <- sd(normalized_slopes, na.rm = TRUE)
# Calculate the coefficient of variation (CV) of slopes
cv_slope <- sd_slope/mean_slope
# Extract p-value for the main effect of the trait
p_value <- summary(model)$coefficients[2, "Pr(>|t|)"]
com_itv <- CV4(subset_data[, paste0(trait, "_minmax_scaled")])
# sp_itv <- mean(tapply(subset_data[, paste0(trait,
# '_minmax_scaled')], list(subset_data$sp), CV4), na.rm = T)
subset_data$new_com <- "Com"
kernel.trait <- TPDs(traits = subset_data[, trait], species = subset_data$new_com) # Kernel density estimates of trait for individuals
functional.indices.trait <- REND(TPDs = kernel.trait)
sp.kernel.trait <- TPDs(traits = subset_data[, trait], species = subset_data$sp)
sp.fun.indices <- REND(TPDs = sp.kernel.trait)
sp_itv <- mean(sp.fun.indices$species$FRichness)
ol <- dissim(sp.kernel.trait)
dis <- ol$populations$dissimilarity
inter_dis <- mean(dis[upper.tri(dis) | lower.tri(dis)])
mean_trait <- mean(subset_data[, paste0(trait, "_minmax_scaled")], na.rm = T)
sp_meantrait <- tapply(subset_data[, paste0(trait, "_minmax_scaled")],
list(subset_data$sp), mean, na.rm = T)
# inter_dis <- CV4(sp_meantrait)
RGR_var <- sd(subset_data$RGR, na.rm = T)
RGR_mean <- mean(subset_data$RGR, na.rm = T)
result <- rbind(result, data.frame(Compete = compete, sample_size = nrow(subset_data),
trait = trait, comp_year = "All", estimate = slope, TLestimate = TLslope,
TMestimate = TMslope, TPestimate = TPslope, TYestimate = TYslope,
Marginal_R2 = r2m, Conditional_R2 = r2c, p_value = p_value, com_ITV = com_itv,
sp_ITV = sp_itv, Fric = functional.indices.trait$species$FRichness,
Fdiv = functional.indices.trait$species$FDivergence, mean_trait = mean_trait,
inter_dis = inter_dis, RGR_var = RGR_var, RGR_mean = RGR_mean, sp_slope_cv = cv_slope))
}
}
return(result)
}
# Core function
CV4 <- function(trait_sample) {
trait_sample <- trait_sample[!is.na(trait_sample)]
if (length(trait_sample) > 1) {
N <- length(trait_sample)
# calcualte CV^2
y_bar <- mean(trait_sample)
s2_hat <- var(trait_sample)
cv_2 <- s2_hat/y_bar^2
cv_1 <- sqrt(cv_2)
gamma_1 <- sum(((trait_sample - y_bar)/s2_hat^0.5)^3)/N
gamma_2 <- sum(((trait_sample - y_bar)/s2_hat^0.5)^4)/N
bias2 <- cv_1^3/N - cv_1/4/N - cv_1^2 * gamma_1/2/N - cv_1 * gamma_2/8/N
cv4 <- cv_1 - bias2
} else {
cv4 <- NA
}
return(cv4)
}
Comp_PCA <- function(alltraits_standard_tda, pr, trait.names, p.title, x.lim, y.lim,
point_shape) {
# Point data for PCA plot
point.pca <- alltraits_standard_tda[, c("Nenv", "sp", "growth_time", "Compete",
"PC1", "PC2")]
# Axis data for PCA plot
x <- "V1"
y <- "V2"
par.all.loading <- data.frame(Trait = trait.names, as.data.frame(matrix(pr$loadings,
nrow = length(trait.names))))
axis.pca <- transform(par.all.loading, v1 = 5 * (get(x)), v2 = 5 * (get(y)))
axis.pca <- data.frame(axis.pca)
# Calculate explained variance for PC1 and PC2
pc1_explained <- round(100 * pr$sd[1]^2/sum(pr$sd^2), 1)
pc2_explained <- round(100 * pr$sd[2]^2/sum(pr$sd^2), 1)
# Create the PCA plot
point.pca1 <- point.pca
point.pca1$cg <- paste0(point.pca1$Compete, point.pca1$growth_time)
point.pca1$cg <- factor(point.pca1$cg, levels = c("alone360", "alone720", "alone1080",
"inter360", "inter720", "inter1080"))
point.pca1$spg <- paste0(point.pca1$sp, point.pca1$growth_time)
point.pca1$spg <- factor(point.pca1$spg, levels = c(unique(point.pca1$spg)))
# point.pca1$compete <- factor(point.pca1$Compete, levels = c('alone',
# 'inter'))
## 绘图
p <- ggplot(point.pca1, aes(x = PC1, y = PC2)) + scale_shape_manual(values = point_shape) +
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ geom_point(aes(color
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ sp,
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ shape
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ Compete),
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ size
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ 1.5,
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ alpha
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ 0.7)
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ +
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ #
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ stat_ellipse(aes(group=sp,color
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ sp),fill=NA,geom='polygon'
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ ,type='norm',alpha=0.2,level
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ 0.95,show.legend
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ T,size=1)+
ggtitle(p.title) + xlab(paste0("PC1 (", pc1_explained, "%)")) + ylab(paste0("PC2 (",
pc2_explained, "%)")) + geom_segment(data = axis.pca, aes(x = 0, y = 0, xend = v1,
yend = v2), arrow = arrow(length = unit(0.2, "cm")), alpha = 0.75, color = "black",
linetype = 2) + geom_text(data = axis.pca, aes(v1 + 0.1, v2 + 0.1, label = Trait),
size = 4, nudge_x = +0.2, nudge_y = +0.3) + theme(legend.text = element_text(face = "italic"),
legend.position = "none") + theme_Publication() + scale_colour_Publication() +
scale_fill_Publication() + xlim(x.lim) + ylim(y.lim)
# Return the PCA plot
return(p)
}
The initial dataset comprised 5188 seedling records, including 1077 dead seedlings. To ensure data quality, individuals with unrealistic growth values were removed. Specifically, seedlings with negative growth due to stem breakage or other causes, as well as those with growth rates exceeding four times the standard deviation, were excluded, resulting in the removal of 41 individuals. After this initial filtering, 4070 seedlings with valid relative growth rate (RGR) data remained.
However, during trait measurement, some seedlings lacked leaf or fine root trait data. To maintain consistency across datasets, these seedlings (a total of 740 individuals) were further excluded. Consequently, the final dataset consisted of 3330 individuals with complete functional trait and growth data. This dataset included 1952 seedlings harvested in the second year and 1378 seedlings harvested in the third year.
# Load processed growth data
load("01_data/processed_data/trait_growth_data.Rdata")
# Filter data for the 2nd and 3rd comp_years and rename RGR column
tda <- tda[tda$comp_year %in% c("2nd", "3rd"), ] %>%
rename(RGR = RGR_H_mean) # 5188 # uses the annual average height growth rate (RGR) as an indicator of plant growth
# Count individuals with memo "S" and remove them
nrow(tda[tda$memo == "S", ]) # 1077
## [1] 1077
tda <- subset(tda, memo != "S") # 4111
# Replace negative growth values with NA and count replacements
sum(sapply(tda["RGR"], function(col) sum(col <= 0, na.rm = TRUE))) # 26
## [1] 26
tda["RGR"] <- lapply(tda["RGR"], function(col) ifelse(col <= 0, NA, col))
# Calculate mean and standard deviation for RGR
mean_rgr <- mean(tda$RGR, na.rm = TRUE)
sd_rgr <- sd(tda$RGR, na.rm = TRUE)
# Identify and count outliers beyond mean ± 4 * SD
sum(tda$RGR > (mean_rgr + 4 * sd_rgr) | tda$RGR < (mean_rgr - 4 * sd_rgr), na.rm = TRUE)
## [1] 3
# Set outlier values to NA
tda$RGR <- ifelse(tda$RGR > (mean_rgr + 4 * sd_rgr) | tda$RGR < (mean_rgr - 4 * sd_rgr), NA, tda$RGR)
# Remove individuals with NA RGR or RGR <= 0
tda <- tda[!is.na(tda$RGR) & tda$RGR > 0, ]
# Extract individuals with all valid traits (non-missing and non-zero values)
tda[trait.names] <- lapply(tda[trait.names], as.numeric) # Convert trait names to numeric
indtraits_tda <- tda %>%
# Select columns
select(sp, sp_abbr, Nenv, ring.num, Light, Moisture, AP, Compete, growth_time, comp_year, BM, RGR, all_of(trait.names)) %>%
mutate(across(
all_of(trait.names[!trait.names %in% c("LTO", "LDMC")]),
~ scale(log10(.))
)) %>%
mutate(across(
all_of(c("LTO", "LDMC")), # "RD" is more normal
~ scale((.))
)) %>%
mutate(across(all_of(trait.names),
~ (. - min(., na.rm = TRUE)) / (max(., na.rm = TRUE) - min(., na.rm = TRUE)),
.names = "{.col}_minmax_scaled"
)) %>%
na.omit() ## 3330
indtraits_tda$RGR <- scale(indtraits_tda$RGR)
# Replace outliers in SLA and LMA with NA
indtraits_tda$SLA[indtraits_tda$SLA < -4] <- NA
indtraits_tda$LMA[indtraits_tda$LMA > 4] <- NA
indtraits_tda$RD[indtraits_tda$RD < -5] <- NA
# Calculate effective sample sizes (non-NA and non-zero values) grouped by year and competition treatment
effective_sample_sizes <- indtraits_tda %>%
group_by(comp_year, Compete) %>%
summarise(across(
all_of(trait.names),
~ sum(!is.na(.) & . != 0),
.names = "{.col}"
)) %>%
pivot_longer(
cols = -c(comp_year, Compete),
names_to = "Trait",
values_to = "Sample_Size"
) %>%
unite("Time_Compete", comp_year, Compete, sep = "_") %>%
pivot_wider(
names_from = Time_Compete,
values_from = Sample_Size
)
kable(effective_sample_sizes)
| Trait | 2nd_alone | 2nd_inter | 3rd_alone | 3rd_inter |
|---|---|---|---|---|
| LA | 1008 | 922 | 692 | 684 |
| SLA | 1008 | 922 | 692 | 681 |
| LMA | 1008 | 922 | 692 | 681 |
| LDMC | 1008 | 922 | 692 | 684 |
| LTO | 1008 | 922 | 692 | 684 |
| SSD | 1008 | 922 | 692 | 684 |
| SRL | 1008 | 922 | 692 | 684 |
| RD | 1008 | 922 | 687 | 680 |
| RTD | 1008 | 922 | 692 | 684 |
| SRA | 1008 | 922 | 692 | 684 |
write.csv(effective_sample_sizes, file = "02_output_results/effective_sample_sizes.csv")
all_indTrait_RGR_results <- data.frame()
for (trait in trait.names) {
trait_result <- calc_trait_RGR(indtraits_tda, trait, "RGR")
all_indTrait_RGR_results <- rbind(all_indTrait_RGR_results, trait_result)
}
indTrait_RGR_results <- data.frame()
for (trait in trait.names) {
trait_annual_result <- calc_annual_trait_RGR(indtraits_tda, trait, "RGR")
indTrait_RGR_results <- rbind(indTrait_RGR_results, trait_annual_result)
}
indTrait_RGR_results$comp_year <- factor(indTrait_RGR_results$comp_year, levels = c("2nd",
"3rd"))
indTrait_RGR_results$trait <- factor(indTrait_RGR_results$trait, levels = trait.names)
all_indTrait_RGR_results$comp_year <- factor(all_indTrait_RGR_results$comp_year,
levels = c("All"))
all_indTrait_RGR_results$trait <- factor(all_indTrait_RGR_results$trait, levels = trait.names)
### Compare competition-free and competition
### Marginal R2
create_ind_trait_plot <- function(indTrait_RGR_results, y, ylab, xlab, title) {
library(ggplot2)
library(ggpubr) # For stat_compare_means
library(patchwork) # If you are using patchwork for layout
library(ggthemes) # For theme_Publication
ggplot(indTrait_RGR_results, aes(x = Compete, y = !!sym(y), group = Compete)) +
geom_boxplot(position = position_dodge(width = 0.6), aes(group = Compete, linetype = Compete), width = 0.6, alpha = 1, size = 0.7) +
geom_jitter(position = position_dodge(width = 0.6), aes(col = trait, group = Compete, shape = Compete), alpha = 0.3, size = 2.5) +
geom_line(aes(group = trait, color = trait), size = 0.8, alpha = 0.3) +
scale_fill_manual(values = c("#cc340c", "#e8490f", "#f18800", "#e4ce00", "#9ec417", "#13a983", "#44c1f0", "#3f60aa", "#739072", "#6F4E37")) +
scale_color_manual(values = c("#cc340c", "#e8490f", "#f18800", "#e4ce00", "#9ec417", "#13a983", "#44c1f0", "#3f60aa", "#739072", "#6F4E37")) +
scale_shape_manual(values = c(1, 19), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
scale_x_discrete(labels = c("alone" = "Competition-free", "inter" = "Competition")) + # 修改 x 轴标签
stat_compare_means(
method = "t.test", paired = TRUE, comparisons = list(c("alone", "inter")),
labels = c("alone" = "Competition-free", "inter" = "Competition")
) +
labs(x = NULL, y = ylab, title = title) +
theme_Publication() + # facet_grid(. ~ comp_year) +
scale_linetype_manual(
values = c("alone" = "dotted", "inter" = "solid"),
labels = c("alone" = "Competition-free", "inter" = "Competition")
) +
theme(legend.position = "none")
}
fig2a <- create_ind_trait_plot(all_indTrait_RGR_results, "Marginal_R2", "Marginal R²", "", "(a)") +
theme_Publication(base_size = 14) +
ylim(c(0, 0.28)) +
theme(legend.position = "none")
all_indTrait_RGR_results$abs_estimate <- abs(all_indTrait_RGR_results$estimate)
fig2b <- create_ind_trait_plot(all_indTrait_RGR_results, "abs_estimate","The absolute value of the slope","","(b)") +
theme_Publication(base_size = 14) +
ylim(c(0, 0.35)) +
theme(legend.position = "none")
combined_fig2_legend <- (fig2a | fig2b) +
plot_layout(nrow = 1, guides = "collect") &
theme(legend.position = "right",
legend.direction = "vertical",
legend.box = "vertical",
legend.spacing.y = unit(0.5, "cm"),
legend.key.size = unit(0.7, "cm"),
axis.line = element_line(size = 1),
axis.ticks = element_line(size = 1),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14))
#combined_fig2_1
### induvidual-level
ind_tda_long <- gather(indtraits_tda, trait, values, trait.names)
ind_tda_long <- ind_tda_long[ind_tda_long$RGR > 0, ]
ind_tda_long <- na.omit(ind_tda_long)
ind_tda_long$comp_year <- factor(ind_tda_long$comp_year, levels = c("2nd", "3rd"))
ind_tda_long$Nenv <- factor(ind_tda_long$Nenv, levels = c("6", "5", "4", "9", "2", "3", "7", "8", "1"))
ind_tda_long$trait <- factor(ind_tda_long$trait, levels = trait.names)
ind_com_fig <- ggplot(ind_tda_long, aes(x = values, y = RGR, col = Compete)) +
geom_point(alpha = 0.1, size = 1, aes(col = Compete)) +
facet_wrap(. ~ trait,
scales = "free", ncol = 5,
labeller = labeller(Compete = c("alone" = "Competition-free", "inter" = "Competition"))) +
geom_smooth(method = "lm", aes(group = Compete, fill = Compete), col = "black", alpha = 0.9, size = 0.15, se = TRUE) +
labs(title = "(c)", x = "Trait values", y = "RGR") +
theme_Publication() +
scale_color_manual(values = c("#6B8A7A", "#DAC0A3"), name = "Treatment", labels = c("alone" = "Competition-free", "inter" = "Competition")) +
scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), name = "Treatment", labels = c("alone" = "Competition-free", "inter" = "Competition")) +
theme_Publication(base_size = 14) +
theme(legend.position = "bottom",
legend.key.size = unit(0.4, "cm"),
legend.text = element_text(size = 15),
legend.title = element_blank()) +
scale_x_continuous(breaks = pretty(ind_tda_long$values, n = 5))
#ind_com_fig
# combined_fig2_1 <- ((fig2a / fig2b) | ind_com_fig) +plot_layout(widths = c(1, 2.5))
combined_fig2ab <- fig2a + fig2b +plot_layout(ncol = 2)
combined_fig2 <- (combined_fig2ab/ ind_com_fig) + plot_layout(heights = c(1.3, 2))
combined_fig2
topptx(combined_fig2_legend, filename = "02_output_results/figs/Fig.2_legend.pptx", height = 4, width = 8)
topptx(combined_fig2ab,filename = "02_output_results/figs/Fig.2ab.pptx", height = 4, width = 8)
ggsave(combined_fig2,filename = "02_output_results/figs/Fig.2_all_R2_slope.png", height = 10, width = 8)
ggsave(combined_fig2,filename = "02_output_results/figs/Fig.2_all_R2_slope.pdf", height = 10, width = 8)
ggsave(combined_fig2,filename = "02_output_results/figs/Fig.2_all_R2_slope.png", height = 10, width = 8)
################
fig_s2a <- create_ind_trait_plot(indTrait_RGR_results, "Marginal_R2", "Marginal R²",
"", "(a)") + theme_Publication(base_size = 14) + facet_grid(. ~ comp_year) +
theme(axis.text.x = element_blank())
indTrait_RGR_results$abs_estimate <- abs(indTrait_RGR_results$estimate)
fig_s2b <- create_ind_trait_plot(indTrait_RGR_results, "abs_estimate", "The absolute value of the slope",
"", "(b)") + theme_Publication(base_size = 14) + facet_grid(. ~ comp_year) +
theme(axis.text.x = element_blank())
combined_fig_s2ab <- fig_s2a + fig_s2b + plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = "bottom", legend.spacing.y = unit(0.15, "cm"), legend.key.size = unit(0.6,
"cm"), axis.line = element_line(size = 1), axis.ticks = element_line(size = 1),
legend.text = element_text(size = 10), legend.title = element_blank())
combined_fig_s2ab
topptx(combined_fig_s2ab, filename = "02_output_results/figs/Fig.S2_R2_slope.pptx",
height = 5, width = 8)
ggsave(combined_fig_s2ab, filename = "02_output_results/figs/Fig.S2_R2_slope.pdf",
height = 4.5, width = 8)
both_indTrait_RGR_results <- rbind(all_indTrait_RGR_results, indTrait_RGR_results)
both_indTrait_RGR_long_re <- both_indTrait_RGR_results %>%
pivot_longer(
cols = c("estimate", "TLestimate", "TMestimate", "TPestimate", "TYestimate"),
names_to = "Type",
values_to = "estimate")%>%
mutate(Type = case_when(
Type == "estimate" ~ "Trait",
Type == "TLestimate" ~ "Trait: Light",
Type == "TMestimate" ~ "Trait: Soil moisture",
Type == "TPestimate" ~ "Trait: Soil phosphorus",
Type == "TYestimate" ~ "Trait: Year"))
both_indTrait_RGR_long_re$abs_estimate <- (abs(both_indTrait_RGR_long_re$estimate))
figS2b <- create_ind_trait_plot(both_indTrait_RGR_long_re[both_indTrait_RGR_long_re$comp_year=="All"&both_indTrait_RGR_long_re$Type!="Trait",], "abs_estimate","The absolute value of the slope","Treament","") + facet_wrap(.~Type,ncol=2,scale="free")+
theme_Publication(base_size = 14) +
theme(legend.position = "bottom",
legend.spacing.y = unit(0.15, "cm"),
legend.key.size = unit(0.5, "cm"),
axis.line = element_line(size = 1),
axis.ticks = element_line(size = 1),
legend.text = element_text(size = 11),
legend.title = element_blank())+guides(
color = guide_legend(title = "trait"), # 显示trait的图例
shape = "none", # 隐藏shape的图例
linetype = "none" # 隐藏linetype的图例
)
figS3 <- create_ind_trait_plot(both_indTrait_RGR_long_re[both_indTrait_RGR_long_re$Type!="Trait", ], "abs_estimate","Absolute estimates in trait—growth models","Treament","") + facet_wrap(.~Type,ncol=2,scale="free")+
theme_Publication(base_size = 14) +facet_grid(comp_year~Type,scale="free") +
theme(legend.position = "bottom",
legend.spacing.y = unit(0.15, "cm"),
legend.key.size = unit(0.5, "cm"),
axis.line = element_line(size = 1),
axis.ticks = element_line(size = 1),
legend.text = element_text(size = 11),
axis.text.x = element_text(angle = 15),
legend.title = element_blank())+guides(
color = guide_legend(title = "trait"), # 显示trait的图例
shape = "none", # 隐藏shape的图例
linetype = "none" # 隐藏linetype的图例
)+ylim(c(0,0.18))
figS3
topptx(figS3, filename = "02_output_results/figs/Fig.S3_interaction_estimate.pptx", height = 9, width = 8)
ggsave(figS3, filename = "02_output_results/figs/Fig.S3_interaction_estimate.pdf", height = 9, width = 8)
topptx(figS2b, filename = "02_output_results/figs/Fig.S2_interaction_estimate.pptx", height = 9.5, width = 8.5)
ggsave(figS2b, filename = "02_output_results/figs/Fig.S2_interaction_estimate.pdf", height = 9.5, width = 8.5)
ggsave(figS2b, filename = "02_output_results/figs/Fig.S2_interaction_estimate.png", height = 8.5, width = 8)
all_indTrait_RGR_results$Significance <- ifelse(all_indTrait_RGR_results$p_value <
0.001, "***", ifelse(all_indTrait_RGR_results$p_value < 0.01, "**", ifelse(all_indTrait_RGR_results$p_value <
0.05, "*", "")))
# Add significance stars based on p_value
indTrait_RGR_results$Significance <- ifelse(indTrait_RGR_results$p_value < 0.001,
"***", ifelse(indTrait_RGR_results$p_value < 0.01, "**", ifelse(indTrait_RGR_results$p_value <
0.05, "*", "")))
indTrait_RGR_results1 <- rbind(all_indTrait_RGR_results, indTrait_RGR_results)
fig_s2a <- ggplot(indTrait_RGR_results1, aes(x = reorder(trait, Marginal_R2), y = Marginal_R2,
fill = Compete)) + geom_bar(stat = "identity", position = "dodge") + facet_grid(. ~
comp_year) + theme_Publication() + coord_flip() + labs(title = "(a)", x = "Trait",
y = expression("Marginal" ~ R^2, "of trait-growth models")) + scale_color_manual(values = c("#6B8A7A",
"#DAC0A3"), name = "Treatment", labels = c(alone = "Competition-free", inter = "Competition")) +
scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), name = "Treatment", labels = c(alone = "Competition-free",
inter = "Competition")) + theme(axis.text.x = element_text(angle = 20, hjust = 1))
fig_s2b <- ggplot(indTrait_RGR_results1, aes(x = reorder(trait, Marginal_R2), y = estimate,
fill = Compete)) + geom_bar(stat = "identity", position = "dodge") + facet_grid(. ~
comp_year) + theme_Publication() + coord_flip() + labs(title = "(b)", x = "",
y = "Estimate of trait in growth models") + scale_color_manual(values = c("#6B8A7A",
"#DAC0A3"), name = "Treatment", labels = c(alone = "Competition-free", inter = "Competition")) +
scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), name = "Treatment", labels = c(alone = "Competition-free",
inter = "Competition")) + theme(axis.text.x = element_text(angle = 20, hjust = 1))
combined_figS3 <- fig_s2a + fig_s2b + plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = "bottom", legend.key.size = unit(0.4, "cm"), legend.text = element_text(size = 12),
legend.title = element_text(size = 12))
combined_figS3
topptx(combined_figS3, filename = "02_output_results/figs/fig_s2_r2_slopes_onlytrait_models.pptx",
height = 5, width = 8)
ggsave(combined_figS3, filename = "02_output_results/figs/fig_s3_r2_slopes_onlytrait_models.png",
height = 5, width = 8)
ind_com_fig + facet_wrap(comp_year ~ trait, scales = "free", ncol = 5, labeller = labeller(Compete = c(alone = "Competition-free",
inter = "Competition")))
fric_indtrait <- ggplot(all_indTrait_RGR_results, aes(x = Compete, y = Fric, group = Compete)) +
geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete), width = 0.5, alpha = 0.5) +
scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
scale_shape_manual(values = c(1, 19), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
stat_compare_means(
method = "t.test", paired = TRUE,
comparisons = list(c("alone", "inter")), labels = c("alone" = "Competition-free", "inter" = "Competition")
) +
labs(x = NULL, y = "Trait variability", title = "(a)") +
theme_Publication() +
scale_x_discrete(labels = c("alone" = "Competition-free", "inter" = "Competition")) +
theme(legend.position = "none")#+ ylim(c(2.8, 5))
rgr_indtrait <- ggplot(indtraits_tda, aes(x = Compete, y = RGR, group = Compete)) +
geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete), width = 0.5, alpha = 0.5) +
scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
scale_shape_manual(values = c(1, 19), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
stat_compare_means(
method = "t.test",
comparisons = list(c("alone", "inter")), labels = c("alone" = "Competition-free", "inter" = "Competition")
) +
labs(x = NULL, y = "Relative growth rate (RGR)", title = "(b)") +
theme_Publication() +
scale_x_discrete(labels = c("alone" = "Competition-free", "inter" = "Competition")) +
theme(legend.position = "none")#+ ylim(c(0, 1.2))
combined.fig3ab <- fric_indtrait + rgr_indtrait + plot_layout(ncol = 2, guides = "collect")&
theme(legend.position = "bottom",
legend.key.size = unit(0.7, "cm"),
legend.text = element_text(size = 13),
legend.title = element_blank())
combined.fig3ab
topptx(combined.fig3ab, filename = "02_output_results/figs/Fig.3_TV_RGR.pptx", height = 4, width = 7)
ggsave(combined.fig3ab, filename = "02_output_results/figs/Fig.3_TV_RGR.pdf", height = 4, width = 7)
itv_indtrait <- ggplot(indTrait_RGR_results, aes(x = comp_year, y = Fric, fill = Compete)) +
geom_boxplot(position = position_dodge(width = 0.7), width = 0.5, alpha = 0.5) +
scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
stat_compare_means(
method = "t.test", paired = TRUE,
comparisons = list(c("alone", "inter")),
label = "p.signif") +
labs(x = "Year under competition", y = "Trait variability", title = "(a)") +
theme_Publication() +
theme(axis.text.x = element_text(angle = 0, hjust = 1),
legend.position = "bottom",
legend.key.size = unit(0.4, "cm"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12)) #+ ylim(c(2.8, 5))
rgr_indtrait1 <- ggplot(indtraits_tda, aes(x = comp_year, y = RGR, fill = Compete)) +
geom_boxplot(position = position_dodge(width = 0.8), width = 0.5, alpha = 0.5) + # 增加 position_dodge 的 width
scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
stat_compare_means(
method = "t.test",
comparisons = list(c("alone", "inter")),
label = "p.signif"
) +
labs(x = "Year under competition", y = "Relative growth rate (RGR)", title = "(b)") +
theme_Publication() +
theme(axis.text.x = element_text(angle = 0, hjust = 1))
combind_figS3_1 <- itv_indtrait + rgr_indtrait1 + plot_layout(ncol = 2, guides = "collect")&
theme(legend.position = "bottom",
legend.key.size = unit(0.7, "cm"),
legend.text = element_text(size = 13),
legend.title = element_blank())
combind_figS3_1
topptx(combind_figS3_1,filename = "02_output_results/figs/Fig.s3_TV_RGR.pptx", height = 4.5, width = 8)
ggsave(combind_figS3_1,filename = "02_output_results/figs/Fig.s3_TV_RGR.pdf", height = 4.5, width = 8)
Fric_indtrait <- ggplot(indTrait_RGR_results, aes(x = Compete, y = com_ITV, group = Compete)) +
geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete),
width = 0.5, alpha = 0.5) + scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"),
labels = c(alone = "Competition-free", inter = "Competition")) + scale_shape_manual(values = c(1,
19), labels = c(alone = "Competition-free", inter = "Competition")) + stat_compare_means(method = "t.test",
paired = TRUE, comparisons = list(c("alone", "inter")), labels = c(alone = "Competition-free",
inter = "Competition")) + labs(x = "Year under competition", y = "Fric",
title = "(a)") + theme_Publication() + facet_grid(. ~ comp_year) + theme(axis.text.x = element_blank()) #+ ylim(c(0.05, 0.6))
Fdiv_indtrait <- ggplot(indTrait_RGR_results, aes(x = Compete, y = Fdiv, group = Compete)) +
geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete),
width = 0.5, alpha = 0.5) + scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"),
labels = c(alone = "Competition-free", inter = "Competition")) + scale_shape_manual(values = c(1,
19), labels = c(alone = "Competition-free", inter = "Competition")) + stat_compare_means(method = "t.test",
paired = TRUE, comparisons = list(c("alone", "inter")), labels = c(alone = "Competition-free",
inter = "Competition")) + labs(x = "Year under competition", y = "FDiv",
title = "(b)") + theme_Publication() + facet_grid(. ~ comp_year) + theme(axis.text.x = element_blank()) #+ ylim(c(0.05, 0.6))
indTrait_RGR_results1 <- rbind(all_indTrait_RGR_results, indTrait_RGR_results)
spitv_indtrait <- ggplot(indTrait_RGR_results1, aes(x = comp_year, y = sp_ITV, fill = Compete)) +
geom_boxplot(position = position_dodge(width = 0.7), width = 0.5, alpha = 0.5) +
scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c(alone = "Competition-free",
inter = "Competition")) + stat_compare_means(method = "t.test", paired = TRUE,
comparisons = list(c("alone", "inter")), labels = c(alone = "Competition-free",
inter = "Competition")) + labs(x = "Year under competition", y = "ITV", title = "(c)") +
theme_Publication() + theme(axis.text.x = element_text(angle = 0, hjust = 1),
legend.position = "bottom", legend.key.size = unit(0.4, "cm"), legend.text = element_text(size = 10),
legend.title = element_text(size = 12))
fig_sp_itv <- ggplot(indTrait_RGR_results, aes(x = Compete, y = sp_ITV, group = Compete)) +
geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete),
width = 0.5, alpha = 0.5) + scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"),
labels = c(alone = "Competition-free", inter = "Competition")) + scale_shape_manual(values = c(1,
19), labels = c(alone = "Competition-free", inter = "Competition")) + stat_compare_means(method = "t.test",
paired = TRUE, comparisons = list(c("alone", "inter")), labels = c(alone = "Competition-free",
inter = "Competition")) + labs(x = "Year under competition", y = "Community trait variability",
title = "(d)") + theme_Publication() + facet_grid(. ~ comp_year) + theme(axis.text.x = element_blank()) #+ ylim(c(0.05, 0.6))
interdis_indtrait <- ggplot(indTrait_RGR_results, aes(x = Compete, y = com_ITV, group = Compete)) +
geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete),
width = 0.5, alpha = 0.5) + scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"),
labels = c(alone = "Competition-free", inter = "Competition")) + scale_shape_manual(values = c(1,
19), labels = c(alone = "Competition-free", inter = "Competition")) + stat_compare_means(method = "t.test",
paired = TRUE, comparisons = list(c("alone", "inter")), labels = c(alone = "Competition-free",
inter = "Competition")) + labs(x = "Year under competition", y = "Community trait variability",
title = "(d)") + theme_Publication() + facet_grid(. ~ comp_year) + theme(axis.text.x = element_blank()) #+ ylim(c(0.05, 0.6))
Fric_indtrait + Fdiv_indtrait + spitv_indtrait + interdis_indtrait + plot_layout(ncol = 2)
# Add the 'growth_year' column based on the 'comp_year' column
indTrait_RGR_results$growth_year <- ifelse(indTrait_RGR_results$comp_year == "2nd", 2,
ifelse(indTrait_RGR_results$comp_year == "3rd", 3, 1))
# View the updated data structure
indTrait_RGR_results$growth_year <- scale(indTrait_RGR_results$growth_year)
ind_trait_model_beta <- glmmTMB(Marginal_R2 ~ com_ITV + com_ITV:Compete + RGR_mean + (1 | trait),
data = all_indTrait_RGR_results,
family = beta_family()
)
summary(ind_trait_model_beta)
## Family: beta ( logit )
## Formula:
## Marginal_R2 ~ com_ITV + com_ITV:Compete + RGR_mean + (1 | trait)
## Data: all_indTrait_RGR_results
##
## AIC BIC logLik deviance df.resid
## -86.5 -80.5 49.2 -98.5 14
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## trait (Intercept) 0.06429 0.2536
## Number of obs: 20, groups: trait, 10
##
## Dispersion parameter for beta family (): 575
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.4400 0.2016 -12.106 < 2e-16 ***
## com_ITV 1.4608 0.7543 1.937 0.0528 .
## RGR_mean 1.0655 0.1555 6.850 7.36e-12 ***
## com_ITV:Competeinter -0.6760 0.5728 -1.180 0.2380
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Load the glmmTMB package for advanced mixed-effects models
library(glmmTMB)
# Fit a Beta regression model using glmmTMB()
ind_trait_model_beta <- glmmTMB(
Marginal_R2 ~ com_ITV + com_ITV:Compete + RGR_mean + growth_year + (1 | trait),
data = indTrait_RGR_results,
family = beta_family()
)
summary(ind_trait_model_beta)
## Family: beta ( logit )
## Formula:
## Marginal_R2 ~ com_ITV + com_ITV:Compete + RGR_mean + growth_year +
## (1 | trait)
## Data: indTrait_RGR_results
##
## AIC BIC logLik deviance df.resid
## -179.0 -167.2 96.5 -193.0 33
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## trait (Intercept) 0.1319 0.3632
## Number of obs: 40, groups: trait, 10
##
## Dispersion parameter for beta family (): 237
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.82905 0.26345 -10.739 < 2e-16 ***
## com_ITV 1.77990 0.98745 1.803 0.0715 .
## RGR_mean 1.13453 0.19170 5.918 3.25e-09 ***
## growth_year 0.22517 0.05346 4.212 2.53e-05 ***
## com_ITV:Competeinter -0.62508 0.68126 -0.918 0.3589
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(ind_trait_model_beta)
## R2m R2c
## [1,] 0.7051071 0.9282304
# Extract coefficients and confidence intervals
coef_ind_cr <- broom.mixed::tidy(ind_trait_model_beta, effects = "fixed", conf.int = TRUE)
coef_ind_cr <- coef_ind_cr[coef_ind_cr$term != "(Intercept)", ]
coef_ind_cr$col <- ifelse(coef_ind_cr$conf.low * coef_ind_cr$conf.high < 0, "white", "black")
colnames(coef_ind_cr)[2] <- "Parameter"
Parameter_names <- c("Trait variability", "Mean growth rate", "Year", "Trait variability:Competition")
# Define parameter levels and plot
coef_ind_cr$Parameter <- factor(Parameter_names, levels = c("Year", "Mean growth rate", "Trait variability:Competition", "Trait variability"))
fig_traitvar <- ggplot(coef_ind_cr, aes(x = Parameter, y = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.4) +
geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3) +
xlab(NULL) +
ylab(expression("Effect size of " * italic("Marginal R"^2) * " in trait-growth")) +
ggtitle("(c)") +
theme_Publication(base_family = "sans") +
coord_flip() +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 3.55, xmax = 4.45), fill = "#6B8A7A", color = NA, alpha = 0.1) +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 2.55, xmax = 3.45), fill = "#DAC0A3", color = NA, alpha = 0.1) +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 0.55, xmax = 2.45), fill = "grey", color = NA, alpha = 0.1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3)
fig_traitvar
# itv_indtrait + rgr_indtrait + fig_traitvar + plot_layout(ncol = 2, nrow = 2, design = "ABCC")
topptx(fig_traitvar,
filename = "02_output_results/figs/Fig3.2_mixed_effect_model.pptx", height = 5, width = 7
)
library(ggrepel)
r_square_itv <- ggplot(indTrait_RGR_results, aes(x = (com_ITV), y = Marginal_R2, color = Compete)) +
geom_point(alpha = 1, size = 2) + # Add points with transparency for better visualization
geom_smooth(method = "lm", size = 1, alpha = 0.5, fill="grey70") +
geom_text_repel(aes(label = trait), size = 3, max.overlaps = 10) +
labs(x = "ITV", y = "Marginal R²", color = "Compete") + # Axis and legend labels
theme_Publication() +
facet_grid(. ~ comp_year) +
stat_poly_eq(
method = "lm",
aes(label = paste(after_stat(rr.label),
after_stat(p.value.label),
sep = "*\", \"*"
)), size = 3
) +
theme(legend.position = "bottom",
legend.key.size = unit(0.4, "cm"),
legend.text = element_text(size = 15),
legend.title = element_blank()) +
scale_color_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition"))
r_square_itv
# Fit a Beta regression model using glmmTMB()
all_ind_trait_model_beta <- glmmTMB(Marginal_R2 ~ Fric + Fric:Compete + RGR_mean + (1 | trait),
data = all_indTrait_RGR_results, family = beta_family())
summary(all_ind_trait_model_beta)
## Family: beta ( logit )
## Formula: Marginal_R2 ~ Fric + Fric:Compete + RGR_mean + (1 | trait)
## Data: all_indTrait_RGR_results
##
## AIC BIC logLik deviance df.resid
## -91.0 -85.0 51.5 -103.0 14
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## trait (Intercept) 0.03922 0.198
## Number of obs: 20, groups: trait, 10
##
## Dispersion parameter for beta family (): 588
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.8241 0.5980 -6.395 1.61e-10 ***
## Fric 0.9575 0.3034 3.156 0.00160 **
## RGR_mean -2.2323 1.2726 -1.754 0.07942 .
## Fric:Competeinter -0.9836 0.3564 -2.760 0.00579 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ind_trait_model_beta <- glmmTMB(
Marginal_R2 ~ Fric + Fric:Compete + RGR_mean +(1 | trait) + (1 | comp_year),
data = indTrait_RGR_results, family = beta_family())
summary(ind_trait_model_beta)
## Family: beta ( logit )
## Formula: Marginal_R2 ~ Fric + Fric:Compete + RGR_mean + (1 | trait) +
## (1 | comp_year)
## Data: indTrait_RGR_results
##
## AIC BIC logLik deviance df.resid
## -175.3 -163.5 94.6 -189.3 33
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## trait (Intercept) 1.203e-01 3.468e-01
## comp_year (Intercept) 5.449e-11 7.382e-06
## Number of obs: 40, groups: trait, 10; comp_year, 2
##
## Dispersion parameter for beta family (): 205
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.00448 0.60718 -4.948 7.49e-07 ***
## Fric 0.29915 0.17226 1.737 0.0825 .
## RGR_mean 0.10839 0.19267 0.563 0.5737
## Fric:Competeinter -0.32091 0.05605 -5.726 1.03e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(ind_trait_model_beta)
## R2m R2c
## [1,] 0.6915189 0.9108413
###################################################
ind_trait_model_beta1 <- glmmTMB(Marginal_R2 ~ Fric + Fric:Compete + (1 | trait),
data = all_indTrait_RGR_results, family = beta_family())
summary(ind_trait_model_beta1)
## Family: beta ( logit )
## Formula: Marginal_R2 ~ Fric + Fric:Compete + (1 | trait)
## Data: all_indTrait_RGR_results
##
## AIC BIC logLik deviance df.resid
## -89.8 -84.8 49.9 -99.8 15
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## trait (Intercept) 0.05413 0.2327
## Number of obs: 20, groups: trait, 10
##
## Dispersion parameter for beta family (): 567
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.29265 0.53643 -6.138 8.35e-10 ***
## Fric 0.49223 0.15401 3.196 0.00139 **
## Fric:Competeinter -0.35972 0.02241 -16.053 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ind_trait_model_beta2 <- glmmTMB(Marginal_R2 ~ RGR_mean + (1 | trait),
data = all_indTrait_RGR_results, family = beta_family())
summary(ind_trait_model_beta2)
## Family: beta ( logit )
## Formula: Marginal_R2 ~ RGR_mean + (1 | trait)
## Data: all_indTrait_RGR_results
##
## AIC BIC logLik deviance df.resid
## -86 -82 47 -94 16
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## trait (Intercept) 0.06781 0.2604
## Number of obs: 20, groups: trait, 10
##
## Dispersion parameter for beta family (): 409
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.17922 0.09136 -23.85 <2e-16 ***
## RGR_mean 1.21163 0.07922 15.29 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Extract coefficients and confidence intervals
coef_ind_cr1 <- broom.mixed::tidy(ind_trait_model_beta1, effects = "fixed", conf.int = TRUE)
coef_ind_cr1 <- coef_ind_cr1[coef_ind_cr1$term != "(Intercept)", ]
coef_ind_cr2 <- broom.mixed::tidy(ind_trait_model_beta2, effects = "fixed", conf.int = TRUE)
coef_ind_cr2 <- coef_ind_cr2[coef_ind_cr2$term != "(Intercept)", ]
coef_ind_cr <- rbind(coef_ind_cr1, coef_ind_cr2)
coef_ind_cr$col <- ifelse(coef_ind_cr$conf.low * coef_ind_cr$conf.high < 0, "white", "black")
colnames(coef_ind_cr)[2] <- "Parameter"
Parameter_names <- c("Trait variability", "Trait variability:Competition","Mean growth rate")
# Define parameter levels and plot
coef_ind_cr$Parameter <- factor(Parameter_names, levels = c("Mean growth rate","Trait variability:Competition", "Trait variability"))
fig3c <- ggplot(coef_ind_cr, aes(x = Parameter, y = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.4) +
geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3) +
xlab(NULL) +
ylab(expression("Effect size of " * italic("Marginal R"^2) * " in trait-growth")) +
ggtitle("(c)") +
theme_Publication(base_family = "sans") +
coord_flip() +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 2.55, xmax = 3.45), fill = "#6B8A7A", color = NA, alpha = 0.1) +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 1.55, xmax = 2.45), fill = "#DAC0A3", color = NA, alpha = 0.1) +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 0.55, xmax = 1.45), fill = "grey60", color = NA, alpha = 0.1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3)
fig3c
combined.fig3 <- (fric_indtrait|rgr_indtrait)/fig3c +
plot_layout(heights = c(1, 0.75))
combined.fig3
topptx(combined.fig3ab, filename = "02_output_results/figs/Fig.3ab_TV_RGR.pptx", height = 4, width = 7)
ggsave(combined.fig3ab, filename = "02_output_results/figs/Fig.3ab_TV_RGR.pdf", height = 4, width = 7)
topptx(fig3c, filename = "02_output_results/figs/Fig.3c.pptx", height = 3, width = 6)
ggsave(fig3c, filename = "02_output_results/figs/Fig.3c.pdf", height = 3, width = 6)
library(ggrepel)
figS3c <- ggplot(both_indTrait_RGR_results, aes(x = Fric, y = Marginal_R2, color = Compete)) +
geom_point(alpha = 1, size = 2) + # Add points with transparency for better visualization
geom_smooth(method = "lm", size = 1, alpha = 0.5, fill="grey70",se=F) +
geom_text_repel(aes(label = trait), size = 3.6, max.overlaps = 10) +
labs(x = "Trait variability", y = "Marginal R²", color = "Compete") + # Axis and legend labels
facet_grid(. ~ comp_year,scale="free") +
stat_poly_eq(
method = "lm",
aes(label = paste(after_stat(rr.label),
after_stat(p.value.label),
sep = "*\", \"*"
)), size = 3
) + theme_Publication(base_size = 14) +
theme(legend.position = "bottom",
legend.key.size = unit(0.4, "cm"),
legend.text = element_text(size = 13),
legend.title = element_blank()) + ylim(c(0,0.34))+
scale_color_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition"))
figS3c
topptx(figS3c, filename = "02_output_results/figs/Fig.S3c.pptx", height = 4, width = 8)
ggsave(figS3c, filename = "02_output_results/figs/Fig.S3c.pdf", height = 4, width = 8)
ind_trait_model_beta <- glmmTMB(Marginal_R2 ~ sp_slope_cv + sp_slope_cv:Compete + (1 | trait),
data = all_indTrait_RGR_results,
family = beta_family()
)
summary(ind_trait_model_beta)
## Family: beta ( logit )
## Formula: Marginal_R2 ~ sp_slope_cv + sp_slope_cv:Compete + (1 | trait)
## Data: all_indTrait_RGR_results
##
## AIC BIC logLik deviance df.resid
## -80.8 -75.8 45.4 -90.8 15
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## trait (Intercept) 0.06259 0.2502
## Number of obs: 20, groups: trait, 10
##
## Dispersion parameter for beta family (): 294
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.2492 0.2641 -8.517 <2e-16 ***
## sp_slope_cv 0.8164 0.3130 2.608 0.0091 **
## sp_slope_cv:Competeinter -1.5147 0.1241 -12.202 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Fit a Beta regression model using glmmTMB()
ind_trait_model_beta <- lmer(
estimate ~ sp_slope_cv + sp_slope_cv:Compete + (1 | trait) + (1 | comp_year),
data = indTrait_RGR_results
)
summary(ind_trait_model_beta)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: estimate ~ sp_slope_cv + sp_slope_cv:Compete + (1 | trait) +
## (1 | comp_year)
## Data: indTrait_RGR_results
##
## REML criterion at convergence: -70.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3867 -0.3881 -0.1127 0.2269 1.9702
##
## Random effects:
## Groups Name Variance Std.Dev.
## trait (Intercept) 0.014616 0.12090
## comp_year (Intercept) 0.000000 0.00000
## Residual 0.003786 0.06153
## Number of obs: 40, groups: trait, 10; comp_year, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.08023 0.06320 32.39836 -1.269 0.21335
## sp_slope_cv 0.19409 0.06682 29.35917 2.905 0.00692 **
## sp_slope_cv:Competeinter -0.06418 0.02484 28.14070 -2.583 0.01527 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sp_sl_
## sp_slope_cv -0.767
## sp_slp_cv:C 0.142 -0.366
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
r.squaredGLMM(ind_trait_model_beta)
## R2m R2c
## [1,] 0.06783866 0.8082099
# Extract coefficients and confidence intervals
coef_ind_cr <- broom.mixed::tidy(ind_trait_model_beta, effects = "fixed", conf.int = TRUE)
coef_ind_cr <- coef_ind_cr[coef_ind_cr$term != "(Intercept)", ]
coef_ind_cr$col <- ifelse(coef_ind_cr$conf.low * coef_ind_cr$conf.high < 0, "white", "black")
colnames(coef_ind_cr)[2] <- "Parameter"
Parameter_names <- c("sp_slope_cv", "sp_slope_cv:Competition")
# Define parameter levels and plot
coef_ind_cr$Parameter <- factor(Parameter_names, levels = c("sp_slope_cv:Competition", "sp_slope_cv"))
fig_traitvar <- ggplot(coef_ind_cr, aes(x = Parameter, y = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.4) +
geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3) +
xlab(NULL) +
ylab(expression("Effect size of " * italic("Marginal R"^2) * " in trait-growth")) +
ggtitle("(c)") +
theme_Publication(base_family = "sans") +
coord_flip() +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 1.55, xmax = 2.45), fill = "#6B8A7A", color = NA, alpha = 0.1) +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 0.55, xmax = 1.45), fill = "#DAC0A3", color = NA, alpha = 0.1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3)
fig_traitvar
library(ggrepel)
r_square_itv <- ggplot(indTrait_RGR_results, aes(x = sp_slope_cv, y = estimate, color = Compete)) +
geom_point(alpha = 0.7, size = 2) + # Add points with transparency for better visualization
geom_smooth(method = "lm", size = 1, alpha = 0.5) +
geom_text_repel(aes(label = trait), size = 3, max.overlaps = 10) +
labs(x = "Estimate", y = "Marginal R²", color = "Compete") + # Axis and legend labels
theme_Publication() +
facet_grid(. ~ comp_year) +
stat_poly_eq(
method = "lm",
aes(label = paste(after_stat(rr.label),
after_stat(p.value.label),
sep = "*\", \"*"
)), size = 3
) +
scale_color_manual(values = c("#609966", "#FFBF61"), labels = c("alone" = "Competition-free", "inter" = "Competition"))
r_square_itv
# Fit a Beta regression model using glmmTMB()
ind_trait_model_beta <- glmmTMB(
Marginal_R2 ~ Fdiv + Fdiv:Compete + (1 | trait) + (1 | comp_year),
data = indTrait_RGR_results,
family = beta_family()
)
summary(ind_trait_model_beta)
## Family: beta ( logit )
## Formula:
## Marginal_R2 ~ Fdiv + Fdiv:Compete + (1 | trait) + (1 | comp_year)
## Data: indTrait_RGR_results
##
## AIC BIC logLik deviance df.resid
## -179.8 -169.7 95.9 -191.8 34
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## trait (Intercept) 7.603e-02 2.757e-01
## comp_year (Intercept) 5.473e-11 7.398e-06
## Number of obs: 40, groups: trait, 10; comp_year, 2
##
## Dispersion parameter for beta family (): 196
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.8551 0.8240 -4.678 2.89e-06 ***
## Fdiv 3.5394 1.5105 2.343 0.0191 *
## Fdiv:Competeinter -2.1671 0.1722 -12.587 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(ind_trait_model_beta)
## R2m R2c
## [1,] 0.7587212 0.9030558
# Extract coefficients and confidence intervals
coef_ind_cr <- broom.mixed::tidy(ind_trait_model_beta, effects = "fixed", conf.int = TRUE)
coef_ind_cr <- coef_ind_cr[coef_ind_cr$term != "(Intercept)", ]
coef_ind_cr$col <- ifelse(coef_ind_cr$conf.low * coef_ind_cr$conf.high < 0, "white", "black")
colnames(coef_ind_cr)[2] <- "Parameter"
Parameter_names <- c("Fdiv", "Fdiv:Competition")
# Define parameter levels and plot
coef_ind_cr$Parameter <- factor(Parameter_names, levels = c("Fdiv:Competition", "Fdiv"))
fig_traitvar <- ggplot(coef_ind_cr, aes(x = Parameter, y = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.4) +
geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3) +
xlab(NULL) +
ylab(expression("Effect size of " * italic("Marginal R"^2) * " in trait-growth")) +
ggtitle("(c)") +
theme_Publication(base_family = "sans") +
coord_flip() +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 1.55, xmax = 2.45), fill = "#6B8A7A", color = NA, alpha = 0.1) +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 0.55, xmax = 1.45), fill = "#DAC0A3", color = NA, alpha = 0.1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3)
fig_traitvar
library(ggrepel)
r_square_itv <- ggplot(indTrait_RGR_results, aes(x = Fdiv, y = Marginal_R2, color = Compete)) +
geom_point(alpha = 0.7, size = 2) + # Add points with transparency for better visualization
geom_smooth(method = "lm", size = 1, alpha = 0.5) +
geom_text_repel(aes(label = trait), size = 3, max.overlaps = 10) +
labs(x = "Fdiv", y = "Marginal R²", color = "Compete") + # Axis and legend labels
theme_Publication() +
facet_grid(. ~ comp_year) +
stat_poly_eq(
method = "lm",
aes(label = paste(after_stat(rr.label),
after_stat(p.value.label),
sep = "*\", \"*"
)), size = 3
) +
scale_color_manual(values = c("#609966", "#FFBF61"), labels = c("alone" = "Competition-free", "inter" = "Competition"))
r_square_itv
# Load the glmmTMB package for advanced mixed-effects models
library(glmmTMB)
# Fit a Beta regression model using glmmTMB()
ind_trait_model_beta <- glmmTMB(
Marginal_R2 ~ sp_ITV + sp_ITV:Compete + (1 | trait) + (1 | comp_year),
data = indTrait_RGR_results,
family = beta_family()
)
summary(ind_trait_model_beta)
## Family: beta ( logit )
## Formula:
## Marginal_R2 ~ sp_ITV + sp_ITV:Compete + (1 | trait) + (1 | comp_year)
## Data: indTrait_RGR_results
##
## AIC BIC logLik deviance df.resid
## -161.3 -151.1 86.6 -173.3 34
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## trait (Intercept) 5.576e-02 2.361e-01
## comp_year (Intercept) 8.042e-11 8.968e-06
## Number of obs: 40, groups: trait, 10; comp_year, 2
##
## Dispersion parameter for beta family (): 107
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.66890 0.46022 -3.626 0.000288 ***
## sp_ITV -0.10419 0.17661 -0.590 0.555245
## sp_ITV:Competeinter -0.39555 0.04996 -7.918 2.42e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(ind_trait_model_beta)
## R2m R2c
## [1,] 0.7593817 0.8508415
# Extract coefficients and confidence intervals
coef_ind_cr <- broom.mixed::tidy(ind_trait_model_beta, effects = "fixed", conf.int = TRUE)
coef_ind_cr <- coef_ind_cr[coef_ind_cr$term != "(Intercept)", ]
coef_ind_cr$col <- ifelse(coef_ind_cr$conf.low * coef_ind_cr$conf.high < 0, "white", "black")
colnames(coef_ind_cr)[2] <- "Parameter"
Parameter_names <- c("Trait variability", "Trait variability:Competition")
# Define parameter levels and plot
coef_ind_cr$Parameter <- factor(Parameter_names, levels = c("Year", "Mean growth rate", "Trait variability:Competition", "Trait variability"))
fig_traitvar <- ggplot(coef_ind_cr, aes(x = Parameter, y = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.4) +
geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3) +
xlab(NULL) +
ylab(expression("Effect size of " * italic("Marginal R"^2) * " in trait-growth")) +
ggtitle("(c)") +
theme_Publication(base_family = "sans") +
coord_flip() +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 1.55, xmax = 2.45), fill = "#6B8A7A", color = NA, alpha = 0.1) +
geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 0.55, xmax = 1.45), fill = "#DAC0A3", color = NA, alpha = 0.1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3)
fig_traitvar
# itv_indtrait + rgr_indtrait + fig_traitvar + plot_layout(ncol = 2, nrow = 2, design = "ABCC")
topptx(fig_traitvar,
filename = "02_output_results/figs/Fig3.2_mixed_effect_model.pptx", height = 5, width = 7
)
library(ggrepel)
r_square_itv <- ggplot(indTrait_RGR_results, aes(x = (com_ITV), y = Marginal_R2, color = Compete)) +
geom_point(alpha = 0.7, size = 2) + # Add points with transparency for better visualization
geom_smooth(method = "lm", size = 1, alpha = 0.5) +
geom_text_repel(aes(label = trait), size = 3, max.overlaps = 10) +
labs(x = "ITV", y = "Marginal R²", color = "Compete") + # Axis and legend labels
theme_Publication() +
facet_grid(. ~ comp_year) +
stat_poly_eq(
method = "lm",
aes(label = paste(after_stat(rr.label),
after_stat(p.value.label),
sep = "*\", \"*"
)), size = 3
) +
scale_color_manual(values = c("#609966", "#FFBF61"), labels = c("alone" = "Competition-free", "inter" = "Competition"))
r_square_itv
ind_trait_model2 <- lmerTest::lmer(Marginal_R2 ~ RGR_mean + com_ITV + sp_ITV + com_ITV:Compete +
sp_ITV:Compete + (1 | comp_year) + (1 | trait), data = indTrait_RGR_results)
tab_model(ind_trait_model2)
| Marginal_R2 | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.09 | -0.02 – 0.20 | 0.121 |
| RGR mean | 0.13 | 0.05 – 0.22 | 0.004 |
| com ITV | 0.15 | -0.05 – 0.35 | 0.142 |
| sp ITV | -0.02 | -0.05 – 0.02 | 0.288 |
| com ITV × Competeinter | -0.08 | -0.23 – 0.07 | 0.277 |
| sp ITV × Competeinter | 0.02 | -0.00 – 0.05 | 0.084 |
| Random Effects | |||
| σ2 | 0.00 | ||
| τ00 trait | 0.00 | ||
| τ00 comp_year | 0.00 | ||
| ICC | 0.77 | ||
| N comp_year | 2 | ||
| N trait | 10 | ||
| Observations | 40 | ||
| Marginal R2 / Conditional R2 | 0.501 / 0.885 | ||
# Figure
coef_ind_cr <- rbind(effectsize(ind_trait_model2)[2:6, ])
colnames(coef_ind_cr)[2:3] <- c("Std_Coefficient", "CI")
coef_ind_cr$col <- "black"
coef_ind_cr$col[coef_ind_cr$CI_low * coef_ind_cr$CI_high < 0] <- "white"
coef_ind_cr$Parameter <- factor(c("Mean RGR", "Community trait variability", "Species trait variability",
"Community trait variability:Competition", "Species trait variability:Competition"),
levels = c("Mean RGR", "Community trait variability:Competition", "Species trait variability:Competition",
"Community trait variability", "Species trait variability"))
fig_traitvar_RGR <- ggplot(coef_ind_cr, aes(x = Parameter, y = Std_Coefficient)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + geom_errorbar(aes(ymin = CI_low,
ymax = CI_high), width = 0.4) + geom_point(shape = 21, fill = coef_ind_cr$col,
lwd = 4, size = 3) + xlab(NULL) + ylab(expression("Effect size of " * italic("Marginal R"^2) *
" in trait-growth")) + ggtitle("") + theme_Publication(base_family = "sans") +
coord_flip() + geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 3.55, xmax = 5.45),
fill = "#6B8A7A", color = NA, alpha = 0.1) + geom_rect(aes(ymin = -Inf, ymax = Inf,
xmin = 1.55, xmax = 3.45), fill = "#DAC0A3", color = NA, alpha = 0.1) + geom_rect(aes(ymin = -Inf,
ymax = Inf, xmin = 0.55, xmax = 1.45), fill = "grey", color = NA, alpha = 0.1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + geom_errorbar(aes(ymin = CI_low,
ymax = CI_high), width = 0.4) + geom_point(shape = 21, fill = coef_ind_cr$col,
lwd = 4, size = 3)
fig_traitvar_RGR
library(piecewiseSEM)
data.sem <- all_indTrait_RGR_results[, c("Compete", "comp_year", "com_ITV", "RGR_mean",
"Fric", "Marginal_R2", "trait")]
data.sem$Compete <- ifelse(data.sem$Compete == "alone", 0, ifelse(data.sem$Compete ==
"inter", 1, NA))
data.sem$Compete <- (as.numeric(data.sem$Compete))
data.sem.na <- data.sem
data.sem <- na.omit(data.sem)
sem_fit_full2 <- psem(lmer(RGR_mean ~ Compete + (1 | trait), na.action = na.omit,
data = data.sem), lmer(com_ITV ~ Compete + (1 | trait), na.action = na.omit,
data = data.sem), lmer(Marginal_R2 ~ Compete * com_ITV + RGR_mean + (1 | trait),
na.action = na.omit, data = data.sem))
# summary(sem_fit_full2, standardized = TRUE, fit.measures = TRUE, rsquare =
# TRUE) Run summary (sem.sum2 <- summary(sem_fit_full2)) sem.sum2$coefficients
plot(sem_fit_full2, node_attrs = list(shape = "rectangle", color = "white", fillcolor = "white"),
digits = 2)
########################
sem_fit_full2 <- psem(lmer(RGR_mean ~ Compete + (1 | trait), na.action = na.omit,
data = data.sem), lmer(Fric ~ Compete + (1 | trait), na.action = na.omit, data = data.sem),
lmer(Marginal_R2 ~ Compete * Fric + RGR_mean + (1 | trait), na.action = na.omit,
data = data.sem))
# summary(sem_fit_full2, standardized = TRUE, fit.measures = TRUE, rsquare =
# TRUE) Run summary (sem.sum2 <- summary(sem_fit_full2)) sem.sum2$coefficients
plot(sem_fit_full2, node_attrs = list(shape = "rectangle", color = "white", fillcolor = "white"),
digits = 2)
ind_sp_fig <- ggplot(ind_tda_long[ind_tda_long$trait %in% c("LTO", "LA", "SLA", "SSD",
"RTD"), ], aes(x = (values), y = RGR, col = sp)) + geom_point(alpha = 1, size = 0.5,
aes(col = sp)) + facet_grid(Compete ~ trait, labeller = labeller(Compete = c(alone = "Competition-free",
inter = "Competition"))) + stat_ellipse(aes(group = sp, fill = sp), geom = "polygon",
alpha = 0.5, color = NA, level = 0.9) + geom_smooth(method = "lm", se = F, alpha = 1,
size = 1, aes(col = sp)) + labs(title = "", x = "Trait values", y = "Relative growth rate (RGR)") +
stat_poly_eq(method = "lm", aes(label = paste(after_stat(rr.label), after_stat(p.value.label),
sep = "*\", \"*")), size = 3) + theme_Publication() + scale_fill_Publication() +
scale_colour_Publication()
ind_sp_fig